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Wavelets,  Fractals,  and  Radial  Basis  Functions 
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Abstract — Wavelets  and  radial  basis  functions  (RBFs)  lead  to 
two  distinct  ways  of  representing  signals  in  terms  of  shifted  basis 
functions.  RBFs,  unlike  wavelets,  are  nonlocal  and  do  not  involve 
any  scaling,  which  makes  them  applicable  to  nonuniform  grids.  De¬ 
spite  these  fundamental  differences,  we  show  that  the  two  types 
of  representation  are  closely  linked  together  . . .  through  fractals. 
First,  we  identify  and  characterize  the  whole  class  of  self-similar 
radial  basis  functions  that  can  be  localized  to  yield  conventional 
multiresolution  wavelet  bases.  Conversely,  we  prove  that  for  any 
compactly  supported  scaling  function  <p{x),  there  exists  a  one¬ 
sided  central  basis  function  p+  ( x )  that  spans  the  same  multireso¬ 
lution  subspaces.  The  central  property  is  that  the  multiresolution 
bases  are  generated  by  simple  translation  of  p_ |_  without  any  dila¬ 
tion.  We  also  present  an  explicit  time-domain  representation  of  a 
scaling  function  as  a  sum  of  harmonic  splines.  The  leading  term  in 
the  decomposition  corresponds  to  the  fractional  splines:  a  recent, 
continuous-order  generalization  of  the  polynomial  splines. 

Index  Terms — Central  basis  functions,  fractals,  fractional 
splines,  localization,  multiresolution,  radial  basis  functions,  re¬ 
finement  filter,  scaling  functions,  self-similarity,  splines,  tempered 
distributions,  two-scale  difference  equation,  wavelets. 


I.  Introduction 

RADIAL  basis  functions  (RBFs)  constitute  a  powerful  tool 
for  working  with  data  that  is  nonuniformly  sampled.  They 
are  used  extensively  for  scattered  data  interpolation  [1],  Thin- 
plate  splines,  in  particular,  have  some  remarkable  variational 
properties:  They  minimize  a  Laplacian  energy  functional  which 
makes  the  solution  invariant  to  rigid-body  coordinate  transfor¬ 
mations  [2],  They  are  often  used  in  medical  imaging  for  the 
estimation  of  deformation  fields  based  on  the  specification  of 
fiducial  points  (or  landmarks)  [3]— [5] .  Radial  basis  functions  are 
also  frequently  applied  to  neural  networks;  they  have  been  pro¬ 
posed  as  an  efficient  mean  for  establishing  a  multidimensional 
mapping  between  some  input  feature  space  and  some  target  re¬ 
sponse  given  a  collection  of  (noisy)  input-output  pairs  (learning 
by  example)  [6],  [7].  Their  use  has  been  justified  through  regu¬ 
larization  theory  in  connection  with  their  energy  minimization 
properties  [8] — [10]. 

A  radial  basis  function  approximation  in  p-dimensions  (x  6 
Rp)  has  the  generic  form 

/(x)  =  ^2  afcKllx-xfcll)  (B 

ke  z 

where  p(r):  R+  — >  R  is  a  univariate  function  and  where 
||x  —  Xfc  ||  denotes  the  Euclidean  distance  between  the  p-vectors 
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x  and  X/. .  The  basis  functions  in  (1)  depend  only  on  the 
distance  to  their  corresponding  grid  point  X/.  and  are  thus  called 
radial.  The  are  weighting  coefficients  that  are  typically 
determined  by  fitting  the  function  to  some  data.  In  the  classical 
interpolation  problem,  the  function  /  is  determined  such  that 
/(xfc)  =  /; . ,  where  the  /;.  are  some  data  values;  in  this  case, 
there  is  exactly  one  linear  constraint  per  basis  function,  and  the 
corresponding  linear  system  of  equations  is  invertible  under 
relatively  mild  conditions  [11].  The  better-known  examples  of 
radial  basis  functions  are  p(r)  =  r  (linear  or  membrane  spline), 
p(r)  =  r2logr  (thin-plate  spline),  and  p(r)  —  \/c2  +  r2 
(Hardy’s  multiquadric). 

Here,  we  are  interested  in  discussing  the  connection  between 
such  radial  basis  functions  and  wavelets.  In  this  latter  context, 
the  generic  multiresolution  approximation  of  a  function  /(a;)  at 
scale  a  —  2~''  has  the  form 

fi(x)  -  ^2  (2) 

ke  z 

where  the  basis  functions  at  level  i  6  Z  are  dilated  by  the 
factor  a  and  spaced  accordingly.  While  (1)  and  (2)  both  involve 
shifts  of  some  basic  template  function — p(||x||),  and  ip(2‘x), 
respectively — there  are  some  crucial  differences  that  need  to 
be  emphasized.  First,  although  it  is  easy  to  define  a  multires¬ 
olution  hierarchy  of  radial  basis  functions  by  simply  dropping 
or  adding  some  points,  there  is  no  provision  for  scaling  in  (1). 
Second,  the  generating  functions  used  in  both  cases  are  fun¬ 
damentally  different:  <p(x)  is  a  relatively  tame,  well-localized 
function  which  is  at  least  square-integrable;  p(r),  on  the  other 
hand,  typically  extends  over  the  whole  axis  and  is  unbounded  at 
infinity.  Conversely,  p{r)  has  a  convenient  closed  form  expres¬ 
sion,  whereas  <p{x)  usually  has  no  such  formula  (it  is  defined 
implicitly  through  a  infinite  recusion)  [12]— [14],  Third,  the  ra¬ 
dial  basis  function  framework  is  ideally  suited  for  a  nonuniform 
setting,  whereas  conventional  wavelet  theory  is  restricted  to  uni¬ 
form  grids. 

Yet,  there  is  an  interesting  link  between  both  types  of 
representations.  This  was  recognized  early  on  by  mathemati¬ 
cians  working  in  approximation  theory  and  has  been  the  basis 
for  several  interesting  generalizations  of  wavelets  [15].  The 
first  step  in  this  realization  was  the  construction  multivariate 
scaling  functions  (or  pre-wavelets)  that  are  linear  combinations 
of  some  radial  basis  functions  [16],  [17].  Utreras  specified 
the  class  of  refinable  radial  basis  functions  as  those  whose 
Fourier  transform  is  p(oj)  —  |tu|-7  [17];  this  is  a  condition 
that  is  unnecessarily  restrictive,  as  we  will  see  here.  The  next 
step,  while  still  working  with  a  uniform  grid,  was  to  extend  the 
multiresolution  concept  without  the  two-scale  relation.  This  led 
to  the  construction  of  nonstationary  wavelets  where  the  basis 
functions  at  distinct  scales  are  no  longer  dilates  of  each  other 
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[18],  [19].  Finally,  there  have  also  been  extensions  of  wavelets 
to  nonuniform  grids  based  again  on  radial  basis  functions  [20], 
[21],  A  univariate  example  of  special  interest  that  falls  into  this 
last  category  are  the  nonuniform  spline  wavelets  [22],  [23]. 

Our  purpose  in  this  paper  is  to  strengthen  the  connection  be¬ 
tween  radial  basis  functions  and  wavelets  even  further.  We  will 
mainly  concentrate  on  the  standard  wavelet  setting  (approxima¬ 
tion  of  univariate  functions  on  a  uniform  grid)  and  approach  the 
problem  from  its  two  different  sides.  First,  in  Section  III,  we  will 
show  how  to  construct  scaling  functions  (or  wavelets)  starting 
from  some  radial  basis  function  p{r),  extending  some  of  the  ear¬ 
lier  results  of  Buhmann  and  Utreras.  In  particular,  we  will  iden¬ 
tify  a  self-similarity  condition  for  p(r)  that  justifies  the  use  of 
the  word  fractal  in  the  title.  We  will  then  characterize  the  whole 
class  of  these  “fractal”  functions  and  show  how  these  can  be  lo¬ 
calized  to  yield  valid  scaling  functions.  We  will  illustrate  our 
theoretical  results  with  the  construction  of  fractional  B-splines 
and  other  related  functions.  In  Section  IV,  we  will  look  at  the 
converse  side  of  the  problem  and  prove  that  any  standard  mul¬ 
tiresolution  analysis  of  L2(R)  can  be  expressed  in  terms  of  cen¬ 
tral  basis  functions,  which  are  a  slight  extension  of  the  radial 
ones.  In  other  words,  we  will  uncover  the  central  basis  function 
that  lies  hidden  behind  any  scaling  function  or  wavelet.  Finally, 
by  combining  our  results,  we  will  present  a  general  parametriza- 
tion  of  all  compactly  supported  wavelets  that  will  provide  us 
with  a  better  understanding  of  their  fractal  nature. 

The  next  section  is  meant  as  an  introduction.  Rather  than 
starting  with  mathematics  right  away,  we  will  consider  a  con¬ 
crete  example  that  contains  most  of  the  flavor  of  what  will  be 
investigated  in  more  generality  in  the  later  sections  of  the  paper. 

II.  Motivation:  The  Example  of  Linear  Splines 

Here,  we  use  piecewise  linear  functions  to  build  a  multires¬ 
olution  analysis  of  L2(R),  but  we  proceed  in  a  nonstandard 
fashion. 

A.  Splines  and  One-Sided  Power  Functions 

The  simplest  conceivable  linear  spline  is  the  one-sided  ramp 
(or  power)  function 

x ,  x  >  0 

(3) 

0,  elsewhere 

which  has  exactly  one  knot  (or  discontinuity)  at  the  origin  [cf. 
Fig.  1(a)],  By  simply  shifting  and  adding  up  such  building 
blocks,  we  can  generate  functions  of  the  form 

s(x)  -  ^2  ak{x  -  xk)+  (4) 

k 

where  the  «/.  are  some  (arbitrary)  weights  and  the  an  in¬ 
creasing  sequence  of  knots.  Since  s(x)  in  (4)  is  a  sum  of  splines, 
it  is  a  spline  as  well  that  inherits  the  discontinuities  of  its  indi¬ 
vidual  constituants.  Hence,  it  is  a  piecewise  linear  function  with 
knots  at  x-k-  Thus,  by  placing  the  knots  at  the  integers  (x-k  =  k), 
we  can  specify  the  basic  space1  of  uniform  splines  of  degree  1 

V0  =  span keZ{p+(x  -  k)}  (5) 

this  paper,  span k(=z{p(x  —  &)}  stands  for  the  function  space  of  all  the 
linear  combinations  of  integer  shifts  of  p(x)  that  are  locally  square  integrable. 


(b) 


Fig.  1.  Example  of  linear  splines,  (a)  One-sided  ramp  x^ ,  which  exhibits 
a  discontinuity  at  the  origin,  (b)  B-spline  of  degree  1.  The  function  in  (b)  is 
localized,  whereas  the  one  in  (a)  is  not. 

where  we  have  set  p+(x)  —  a:+. 

By  taking  the  second  forward  finite  differences  of  p+(x) — a 
specific  instance  of  (4)  with  ao  =  1,  cq  =  —2,  and  a 2  =  1 — we 
obtain  the  hat  function  (or  causal  B-spline  of  degree  1) 

(3\{x)  —  A2  x+  —  x+  —  2(x  —  1)+  +  ( x  —  2)+  (6) 

which  is  the  more  standard,  compactly  supported  basis  function 
for  the  linear  splines  [cf.  Fig.  1(b)].  The  remarkable  fact  with 
(6)  is  that  the  one-sided  functions  cancel  out  for  x  >  2;  in  other 
words,  we  are  able  to  localize  p. |_  by  taking  a  suitable  linear 
combination  of  its  shifted  versions.  Moreover,  we  can  invert  (6) 
and  express  the  one-sided  power  function  as  a  weighted  sum  of 
B-spline  basis  functions 

x+  —  ^  (k  +  1  )(3\{x  -  k).  (7) 

k>  0 

This  shows  that  our  definition  (5)  of  the  basic  spline  space  is 
equivalent  to  the  standard  one  that  involves  linear  combinations 
of  B-splines.  Thus,  {p+(a:  —  k)}ke z  is  a  valid  basis  for  Vo, 
albeit  not  a  Riesz  basis;  the  main  difficulty  is  that  the  in 
(4)  are  not  necessarily  in  f2,  which  also  means  that  this  basis 
is  not  well  conditioned.  Thus,  we  will  need  to  exert  special  care 
while  manipulating  expansions  such  as  (4)  both  numerically  and 
mathematically. 

B.  How  Multiresolution  Becomes  Trivial 

The  advantage  of  the  present  formulation  is  that  it  makes 
the  multiresolution  structure  of  splines  stand  out  quite  naturally 
(cf.  Fig.  2).  Consider  the  coarse-to-fine  sequence  of  subspaces 
•  •  •  V_i  C  Vo  C  Vi  ■  •  ■  C  V,:  ■  •  where  V,-  represents  the  space 
of  linear  splines  with  knots  at  x^  —  k2~\  k  6  Z.  These  splines 
are  generated  simply  by  dropping  all  the  basis  functions  in  (5) 
that  are  not  positioned  at  the  desired  knots.  Thus,  we  define  our 
uniform  spline  space  with  scale  T,  —  2~'1  as 

Vi  =  span{p+(a:  -  k2~1)}.  (8) 

kei 
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Fig.  2.  Multiresolution  spaces  using  one-sided  power  functions. 


complement  of  V,  in  V,+i.  Such  nonuniform  spline  wavelets 
were  first  described  by  Lyche  as  early  as  in  1992  [22], 

D.  Central  Basis  Functions 

Until  now,  we  have  worked  with  the  one-sided  function 
p+(x)  —  x+.  We  could  have  also  used  its  reversed  version 
(>- (x)  —  p+{—xf,  which  can  be  localized  using  backward 
finite  (or  divided)  differences.  Another  option  would  have 
been  to  use  the  symmetric  function  (>,.  (x)  =  .r|,  which  is 

the  prototypical  example  of  a  radial  basis  function.  In  the 
sequel,  we  will  consistently  use  the  subscripts  *,  +,  and  —  to 
denote  symmetric  (or  radial),  causal,  and  anti-causal  generating 
functions,  respectively.  These  are  all  special  cases  of  what  we 
call  central  basis  functions  and  which  we  denote  by  the  generic 
symbol  p{x).  The  one-sided  functions  are  typically  simpler  to 
work  with  in  the  univariate  case.  The  symmetric  ones,  on  the 
other  hand,  are  usually  preferred  in  a  multivariate  setting,  for 
they  fall  within  the  framework  the  well-developed  theory  of 
radial  basis  functions. 


Clearly,  the  basis  functions  for  V,  are  a  subset  of  those  of  V’,' 
for  i'  >  i,  which  implies  that  V)  C  V,/,  for  all  i!  >  i,  which 
is  a  multiresolution  property.  Since  each  V,  also  has  a  B-spline 
Riesz  basis  {2^2  rl^(2'  x  —  k)}kez,  the  whole  ladder  of  spline 
subspaces  for  i  6  Z  generates  a  multiresolution  of  L2(F?)  in 
the  sense  specified  by  Mallat  [24].  Hence,  it  is  possible  to  con¬ 
struct  a  whole  variety  of  corresponding  wavelet  bases  using  any 
of  the  standard  design  techniques.  Specific  examples  of  linear 
spline  wavelets  that  are  orthogonal  [24],  semi-orthogonal  [25], 
biorthogonal  [26],  or  even  shift-orthogonal  [27]  have  been  de¬ 
scribed  in  the  literature. 

C.  Nonuniform  Linear  Splines 

The  power  of  the  present  formulation  really  becomes 
apparent  if  we  move  one  step  further  and  consider  a  given 
nonuniform  sequence  of  knots  •  •  •  <  <  •  •  •  with 

k  6  Z.  We  then  define  a  corresponding  embedded  sequence  of 
nonuniform  spline  spaces 

V-i  =  span{p+  (x  -  xk2i )}  (9) 

kez 

which  share  the  same  inclusion  properties  as  before:  •  •  •  C 
V-i  C  . . .  V_i  C  Vo-  Here,  too,  we  are  able  to  produce  com¬ 
pactly  supported  basis  functions  (nonuniform  B-splines),  except 
that  they  lose  the  convenient  shift-invariant  structure  that  is  in¬ 
herent  to  standard  (uniform)  multiresolution  analysis.  Specifi¬ 
cally,  the  fcth  (nonuniform)  B-spline  at  resolution  —i  is 

/  ( X  -  Xk2i)+  ~  -  X(k+1)2i)  +  \ 

y  x(k+l)2i  ~  xk2i  J 

(X  ~  X(k+l)2i)+  ~  (*  ~  ^+2)2’-)  A 

a’(fc+2)2!'  -  x(k+l)2i  J 

It  is  a  triangular  function  that  takes  the  value  one  at  x  —  c  /,■ + 1  ;  2 ■ 
and  vanishes  for  x  <  xk2i  and  x  >  X(k+2-)2i .  Note  that  this  lo¬ 
calization  process  involves  divided  differences  of  p+(x)  rather 
than  finite  differences,  as  in  (6).  This  nonuniform  setting  is  also 
suitable  for  constructing  wavelet  bases  that  span  the  orthogonal 


E.  Extension  to  Higher  Dimensions:  Variational  Splines 

Let  us  now  briefly  show  how  to  extend  this  type  of  construc¬ 
tion  to  a  nonuniform  grid  in  p  dimensions.  Here,  it  becomes 
preferable  to  work  with  the  radial  basis  function  (>+  (x)  —  .r|. 
We  also  need  to  specify  our  nonuniform  grid  at  the  finest  scale 
available.  It  is  given  in  the  form  of  an  ordered  set  of  (distinct) 
points  6  F?p,  with  k  6  Z.  These  may,  for  example,  corre¬ 
spond  to  the  discrete  locations  at  which  the  samples  of  an  input 
function  are  given.  We  then  specify  our  basic  radial  spline  space 

V0  =  span{p*(||x  -  Xfc||)}.  (11) 

kez 

For  [}  >  1,  the  functions  in  Vo  are  no  longer  piecewise  linear. 
Yet,  they  are  compatible  with  the  variational  formulation  of 
linear  splines.  The  problem  is  the  following:  Given  the  samples 
of  a  function  /(x;.)  =  // . ,  k  6  Z,  one  wishes  to  find  the  best 
interpolant  in  L2(F?P)  such  that  an  adequate  Duchon  semi-norm 
is  minimized  [2],  It  can  be  shown  that  the  optimum  solution  is 
necessarily  included  in  Vo;  specifically,  it  is  the  unique  function 
s(x)  £  Vo  that  satisfies  the  interpolation  constraint  s(xk)  =  fk. 
k  £  Z.  The  functional  f  ||  V/(x)||2  dx  may  be  thought  of  as 
the  bending  energy  of  a  thin  membrane;  hence,  we  have  the 
term  membrane  spline.  Finding  the  expansion  coefficients  for 
the  basis  p*(||x  —  X/. ||)  is  precisely  what  the  radial  basis  func¬ 
tions  interpolation  problem  is  all  about. 

F.  Road  Map  to  the  Paper 

The  manipulations  that  we  have  made  so  far  using  linear 
splines  are  also  valid  for  higher  order  splines.  This  is  all  well 
known  in  approximation  theory.  In  the  remainder  of  the  paper, 
we  will  generalize  these  ideas  further  and  show  that  the  first 
part  of  the  process  described  above  (Sections  II-A  and  II-B,  in 
particular)  is  applicable  in  a  much  wider  setting  as  had  been  re¬ 
alized  so  far.  First,  we  will  identify  the  whole  class  of  central 
basis  functions  p(x)  (one-sided  or  radial)  that  can  be  used  to 
specify  a  multiresolution  analysis  of  L2(IR)  in  the  sense  defined 
by  Mallat.  We  will  see  that  an  admissible  central  basis  func¬ 
tion  must  satisfy  a  self-similarity  relation.  This,  together  with 
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a  suitable  localization  process  [similar  to  (6)],  will  allow  us  to 
specify  new  classes  of  scaling  functions  much  more  explicitly 
than  is  done  usually  by  means  of  an  infinite  recursion  ftwo-scale 
relation).  Second,  we  will  look  at  the  converse  problem  and 
prove  that  a  central  basis  function  p(x)  can  be  found  for  any 
compactly  supported — and  more  generally,  one-sided — scaling 
function  ip(x).  The  next  logical  step,  which  will  not  be  pursued 
here,  is  to  use  these  results  to  extend  the  standard  wavelet  con¬ 
structions  to  the  nonuniform  case  similar  to  what  we  have  done 
with  the  linear  splines  in  the  example  above. 

III.  From  Central  Basis  Functions  to  Wavelets 

We  will  now  show  how  it  is  possible  to  construct  wavelets 
starting  from  a  radial  (or  central)  basis  function  p(x).  We 
will  not  consider  wavelets  literally  but,  rather,  their  associated 
scaling  functions,  which  are  the  key  to  the  multiresolution 
structure  of  the  wavelet  transform.  Once  the  scaling  function 
has  been  specified,  it  is  comparatively  easy  to  construct  a 
corresponding  wavelet  basis  using  standard  techniques  [12], 

[13]. 

A.  Notation 

We  will  need  to  raise  complex  numbers  to  some  complex 
power.  Euler  showed  that  this  is  a  well-defined  operation,  as 
long  as  one  specifies  the  argument  in  a  unique,  unambiguous 
fashion,  e.g.,  z  —  |^|eJ  arg(*)  with  arg(z)  6  [— 7t,  7t[.  Then,  we 
specify 

z*  —  e*'(log  1*1+3  argi). 

For  instance,  we  have  \zz' \  =  \z\R’^zJe~'-T^zJax%z . 

Euler’s  Gamma  function  is  defined  for  x  >  0  by 

poo 

r(:r)  =  / 

Jo 

This  is  a  C°°  interpolation  of  the  factorial  since  n!  =  T(n  +  1) 
if  n  6  N.  It  is  further  possible  to  extend  this  function  by  ana¬ 
lytic  continuation  to  complex  values  of  x,  except  at  the  negative 
integers,  where  T(x)  is  unbounded  (see  [28]). 

We  denote  by  f(to)  —  fR  f(x)e~3ZJX  dx  the  Fourier  trans¬ 
form  of  the  absolutely  integrable  function  f(x).  For  functions 
such  as  p(x)  that  are  not  bounded  at  oo  but  do  not  increase  faster 
than  a  polynomial  [i.e.,  p(x)  £  S'  where  S'  is  Schwartz’  class  of 
tempered  distributions],  we  will  use  a  weaker  definition.  Specif¬ 
ically,  p(uj)  is  the  Fourier  transform  of  p(x)  in  the  sense  of  dis¬ 
tributions  if  and  only  if  (ji.  m)  —  (p,  F)  for  any  indefinitely 
differentiable  fastly  decreasing  function  x!  £  S  (“test  function”) 
[29], 

B.  Scaling  Functions 

Often,  a  scaling  function  is  defined  indirectly  through  its  re¬ 
finement  filter  G  [cf.  (12)  below].  One  then  has  to  worry  about 
the  delicate  issues  of  the  convergence  of  the  iterated  filterbank 
and  the  L2(  R)-completcncss  of  the  wavelet  expansion.  Here,  we 
propose  a  more  explicit  definition  that  avoids  these  problems  at 
the  outset. 

Definition  1:  p(x)  is  a  valid  scaling  function  if  and  only  if  it 
satifies  the  following  three  conditions. 


l)  Two-Scale  Relation:  The  function  <p  at  a  coarser  resolu¬ 
tion  can  be  expressed  using  its  shifts 

<p{x/2)  -  22  gk<p{x  -  k).  (12) 

kez 

ii)  L2(R)  Stability:  {ip(x  —  k)}ke z  is  a  Riesz  basis  of  Vo, 
which  is  equivalent  to  require  that  there  exist  two  positive 
constants  A  and  B  such  that 

A  <  \Jp{fi>  +  27t/;)|2  <  B  (13) 

kez 

almost  everywhere. 

iii)  Partition  of  Unity: 

'22  P(x  —  k)  —  1  <£(27 rk)  —  8 (14) 

kez 

The  two  first  conditions  are  implicit  to  Mallat’s  axiomatic 
definition  of  a  multiresolution  analysis  of  L2(R)  [24],  More¬ 
over,  it  can  be  shown  that  Mallat’s  completeness  requirement 
[denseness  in  L2(R)]  is  equivalent  the  partition  of  unity  condi¬ 
tion,  which  has  the  merit  of  being  explicit  (cf.  [30]).  It  follows 
that  the  three  conditions  are  sufficient  to  build  a  multiresolution 
analysis  of  L2(R). 

C.  Admissible  Central  and  Radial  Basis  Functions 

A  central  (either  one-sided  or  radial)  basis  function  p(x)  has 
properties  that  are  quite  distinct  from  those  of  a  scaling  function 
tp  because  of  fundamental  differences  in  the  replication  mech¬ 
anism.  With  radial  basis  functions,  the  multiresolution  bases 
(e.g.,  { p(x  —  27>:)};.cz)  are  generated  by  simple  translation  of 
p{x),  as  opposed  to  dilation  plus  translation  (e.g.,  {p(x/21  — 
k)}ke z)>  as  is  traditionally  the  case  with  the  wavelet  transform. 
The  most  evident  difference  between  the  two  species  of  func¬ 
tions  is  that  the  central  basis  functions  are  not  well  localized.  To 
avoid  gaps,  they  must  be  infinitely  supported.  In  general,  they 
are  not  even  integrable,  and  yet,  it  can  be  shown  that  this  ap¬ 
parent  disadvantage  brings  approximation  order  (i.e.,  the  ability 
to  reproduce  polynomials)  to  the  generated  spaces. 

Definition  2:  A  central  basis  function  p(x)  £  S'  is  admis¬ 
sible  if  and  only  if  the  series  of  subspaces  V,:  =  span/i,sZ{p(a<  — 
k2~')}  for  /  £  Z  are  such  that  we  have  the  following. 

i)  Vi  is  localizable,  i.e.,  V,-  IT  L2(R)  {0}. 

ii)  There  exists  :.p,  £  V,  FI  L2(R)  such  that 

—  L2(R)  stability:  {ipfix  —  k2~t)}kez  is  a  Riesz  basis 

of  v*  n  l2(r) 

—  Partition  of  unity:  pfix  —  k2~‘  )  —  1. 

The  existence  of  i pi  £  V,  FI  L2(R)  means  that  p(x)  can  be 
localized  by  taking  a  suitable  linear  combination  of  its  shifts. 
The  implication  of  the  partition  of  unity  is  more  intriguing.  This 
means  that  it  is  possible  to  represent  the  constant  as  a  linear 
combination  of  p(x  —  kTi ),  irrespective  of  the  scale  l  t  =  2~‘ . 
This  nonunique  way  of  writing  the  constant  implies  that  p{x) 
itself  cannot  generate  a  Riesz  basis,  which  points  out  another 
fundamental  difference  between  p(x)  and  <p(x). 

By  construction,  the  spaces  V,  are  embedded 

•  •  •  C  V_1  C  Vo  C  •  •  •  Vi  C  Vi+1  C  •  •  • 
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but  this  does  not  mean  that  the  p,  in  Definition  2  are  linked  to 
each  other  through  a  two-scale  relation. 

Some  good  examples  of  admissible  central  basis  functions 
are  the  one-sided  power  functions  xf/n\  with  n  integer.  These 
can  be  localized  using  the  backward  (n  +  l)-order  difference 
operator  to  yield  the  causal  B-splines  of  degree  n  (cf.  [31],  [32]): 

/??(*)  =  A£+1^±  M"1.  (15) 

n\  \  jlo  J 

The  B-splines  satisfy  the  three  conditions  in  Definition  1;  they 
constitute  one  of  the  earliest  example  of  valid  scaling  functions 
[24], 

Localizing  less  simple  central  basis  functions  turns  out  to 
be  more  intricate.  It  is  theoretically  possible  in  the  good  cases 
(see  Buhmann  in  [33])  to  localize  them  by  building  the  corre¬ 
sponding  interpolator,  although  it  might  not  be  an  easy  task  to 
get  the  coefficients  of  the  interpolating  filter,  especially  when 
the  Fourier  transform  of  p  is  not  a  true  function.  Moreover,  the 
interpolator  is  not  always  stable;  see  the  classical  example  of 
the  quadratic  splines  for  which  (> +(->:)  —  x\.  We  will  give  our 
alternative  localization  technique  in  Section  III-F. 


D.  Self-Similar  Central  Basis  Functions 

We  will  say  that  a  function  p  is  self-similar  (for  dyadic 
zooming  factors)  iff  it  satisfies 

p(x)  =  A p  (16) 

for  all  x  6  R. 

Theorem  1:  Every  self-similar  admissible  central  basis  func¬ 
tion  generates  a  valid  multiresolution  of  L2(R). 

In  Section  IV,  we  will  see  that  the  converse  is  also  true, 
provided  that  the  scaling  function  p  is  one-sided  or  compactly 
supported. 

Proof:  Let  p  be  an  admissible  central  basis  function  satis¬ 
fying  (16).  Then,  there  exists  p  6  Vo  which  is  a  Riesz  basis  for 
Vo  IT  L2(R).  Since  p  is  in  the  span  of  p(x  —  k),  p(x/2)  is  in  the 
span  of  p{x/ 2  —  k),  that  is,  in  the  span  of  p(x.  —  2 n)  because  of 
the  self-similarity  of  p.  This  also  implies  that  p(x/2)  is  in  the 
span  of  p{x  —  n).  Thus,  p{x/2)  6  Vo  0  L2(R).  Since  p  is  a 
basis  of  Vo  (T  L2(R),  p{x/2)  is  itself  in  the  span  of  p(x  —  k), 
which  proves  that  p  satisfies  the  scale  relation  (12).  All  the  other 
required  properties  are  provided  by  the  definition  2  of  an  admis¬ 
sible  central  basis  function.  ■ 

Because  of  the  self-similarity  relation,  the  graph  of  such  cen¬ 
tral  basis  functions  are  fractals. 

Remarks: 

•  Equation  (16)  implies  that  p{x)  is  either  0  or  tends  to  in¬ 
finity  at  one  end  ( x  =  0  or  x  — >  oo).  If  |  A|  >  1,  then  p(x) 
tends  to  grow  globally  from  0  to  +oq;  otherwise,  when 
|A|  <  1,  it  globally  tends  to  decrease  from  +<do  to  0. 

•  In  order  for  p  to  be  square  integrable  in  [0,  1],  it  is  neces¬ 
sary  that  |A|  >  l/s/2.  This  is  because 


[  p(x)2  dx  —  —  I  p{2x)2  dx  —  — j  I  p(x)2  dx 

Jo  A-  J o  2A-  J o 

l  r1  ^ 

>  2A2  /o  dX 


•  A  self-similar  central  basis  functions  is  entirely  and 
uniquely  specified  by  its  definition  within  the  interval 
[1/2,  1]  and  its  eigenvalue  A.  More  precisely,  if  we  denote 
by  P[i/2,  i]  the  restriction  of  p  to  [1/2,  1],  we  have 

P&)  =  ^,  (J)- 

ie  i 

Conversely,  this  expression  defines  a  self-similar  central 
basis  function  with  eigenvalue  A. 

•  The  product  and  the  convolution  of  two  self-similar  func¬ 
tions  are  self-similar  as  well.  Admissibility,  however,  does 
not  carry  over  so  easily,  e.g.,  the  product  of  the  admissible 
central  basis  functions  |a:|  x  |a:|  yields  x2,  which  is  not 
admissible. 

E.  General  Parametrization 

A  nice  feature  of  one-sided  self-similar  central  basis  func¬ 
tions  is  that  they  can  be  expressed  as  a  (countable)  sum  of  one¬ 
sided  power  functions. 

Theorem  2:  Consider  a  self-similar  one-sided  central  basis 
function  p+  with  eigenvalue  A.  Then,  p  is  square  integrable  in 
[1,  2]  if  and  only  if  it  can  be  expressed  as 

P+(x)  =  £  7„.^A+J(2"r/k*2)  a.e.  (17) 
nez 

where 

in  =  ^2  /"  ?+(or log2  A_ w(2,!7r/  log 2)  de  as) 

is  a  square  summable  sequence. 

Likewise,  under  the  same  hypothesis,  a  self-similar  sym¬ 
metric  central  basis  function  p*  with  eigenvalue  A  can  be 
expressed  as  a  sum  of  symmetric  monomials 

p*(x)  =  Y,  7»kfog2A+i(2,!7r/log2)  a.e.  (19) 
nez 

where  {7„}„sz  €  ?2  ■ 

For  completeness,  we  have  derived  in  Appendix  A  the  Fourier 
transform  of  (17),  which  is  defined  in  the  sense  of  distributions. 

Proof:  As  we  have  seen,  a  one-sided  dyadic  central  basis 
function  is  completely  defined  by  its  restriction  to  the  interval 
[!>  2]- 

We  define  po(x)  —  x  (los V los2ip+(^x)  and  observe  that 
Po(2x)  is  a  1-periodic  function,  i.e.,  po(2*+1)  =  po(2*)- 
Thus,  we  can  apply  Fourier’s  theorem  to  po(2x).  If  po(x) 
is  square  integrable  in  [1,  2]  (which  is  the  case  when 
p(x)  is  square  integrable  in  [1,  2]),  so  is  pq(2x),  and 
then,  po(2*)  =  J2neZ  7 ne->2n7VX  almost  everywhere,  with 
7„  =  Jq1  po(2^)e_j2,!71'^  d£.  A  change  of  variables  then  yields 
the  expansion  formula  (17). 

When  the  central  basis  function  is  symmetric,  it  can  be  ex¬ 
pressed  as  p*(x)  —  p+(x)  +  p+(—x),  where  p+  is  the  restric¬ 
tion  of  p  to  [0,  +oo[.  Applying  (17)  to  p+  and  noticing  that 
|a:|“  =  +  (— x)f  immediately  yields  (19).  ■ 

Of  course,  we  may  also  consider  more  general  self-similar 
central  basis  functions  and  express  them  as  a  sum  of  two 
self-similar  one-sided  central  basis  functions  pi  and  (>■> : 
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p{x)  —  piix)  +  p2(—x),  each  of  which  can  be  decomposed 
according  to  (17). 

F.  Localization 

Theorem  5:  Let  p+  be  a  one-sided  self-similar  central  basis 
function  with  eigenvalue  A,  |A|  >  l/\/2.  If  its  representation  in 
(17)  is  a  finite  sum,  then  the  filter  P  defined  through  its  z- trans¬ 
form  P(z)  -  J2ke z  Pt-Z~k  as 

pOO 

P(z )-'=  /  p+iQe-V-*-1* 

JO 

poo  tn 

=  £*""/  (20) 

n>  0  n' 

localizes  the  central  basis  function  p+,  i.e.,  <p{x)  — 
PkP+(x  —  k )  is  in  L2(IR). 

Similarly,  if  p  is  a  symmetric  self-similar  central  basis  func¬ 
tion  and  if  its  representation  in  (19)  is  a  finite  sum,  then  the 
symmetric  filter  P  defined  through  its  ^-transform  as 


[34].  This  also  means  that  the  scaling  function  will  reproduce  all 
polynomials  of  degree  N  <  ["log?  A]  and  that  the  corresponding 
wavelet  will  have  [log?  A]  +  1  vanishing  moments  [35], 
Contrary  to  its  inverse,  the  filter  P(z)  in  Theorem  3  is  not  nec¬ 
essarily  causal;  this  means  that  the  localized  function  might  not 
be  one-sided.  However,  when  Piz)  has  only  a  finite  number  of 
poles  {Ck}k=i-..N  outside  the  unit  circle,  it  is  possible  to  choose 
po(z)  =  p(z)  nf=i(i  —  1)  as  an  alternative  localization 

filter;  the  advantage  is  that  Pq  is  causal  and,  thus,  that  ■■fix')  is 
one  sided. 


G.  Refinement  Filter  and  Perfect  Reconstruction  Filterbank 

Once  the  localization  filter  Pie3 UJ)  has  been  specified,  we 
may  obtain  the  Fourier  transform  of  refinement  filter  g  in  (12) 
by  evaluating  the  ratio  between  2<p{2 u>)  and  With  the  hy¬ 

potheses  of  Theorem  3  [finite  sum  representation  in  (17)],  this 
simplifies  to 


G{e3W) 


1  P{e32w) 

A  P(e-A')  ' 


(22) 


/OO 

/a*(£)  ^  j  d£, 

poo  pn 

=  £>""  +  *")/  (21) 

n>0  U- 

localizes  the  symmetric  central  basis  function  p* . 

The  proof  of  this  result  is  given  in  Appendix  B.  The  nonlocal 
character  of  p(x)  is  entirely  due  to  the  O(l/u21+log2  A)  singu¬ 
larity  of  its  Fourier  transform  at  the  origin  (cf.  Appendix  A). 
Thus,  the  essence  of  the  localization  process  is  to  cancel  out 
these  singularities.  The  intuitive  explanation  of  Theorem  3  is 
as  follows.  Since  1  —  e~3M  ~  ju>  as  u>  — i  0,  the  integral  on 
the  first  line  of  (20)  for  z  —  e3M  and  u>  close  to  zero  is  essen¬ 
tially  equivalent  to  the  Fourier  transform  of  p{x).  This  means 
that  the  localization  filter  P(e-,UJ)  will  have  multiple  zeros  at 
oj  —  0,  which  will  precisely  cancel  all  the  singularities  of  p{u>). 
An  ideal  localization  [i.e.,  ip{x)  -a  (/(a:)]  would  be  achieved  if 
we  could  design  a  filter  such  that  <p{u>)  —  p(uj)P{e3IJJ)  ~  1  over 
the  whole  frequency  axis.  While  we  can  obviously  not  enforce 
this  constraint  everywhere  because  P(e 3^)  is  27t -periodic  and 
p(u>)  is  not,  we  make  sure  that  it  is  satisfied  around  the  origin. 
In  effect,  our  localization  method  is  akin  to  using  the  bilinear 
transform  to  convert  a  continuous-time  filter  into  a  discrete  one. 

We  conjecture  that  Theorem  3  still  holds  true  when  (17)  is 
not  finite,  at  least  under  wide  conditions  on  p+(x).  Our  proof, 
however,  can  hardly  be  extended  since,  in  that  case,  p+(u>)  is  no 
longer  a  function,  which  makes  the  product  with  P(e3^')  highly 
questionable. 

Note  that  by  construction,  ip(oj)  is  boundedly  differentiable 
at  0  and  has  singularities  of  order  1  +  log?  A  at  w  —  2mv  for 
n  6  Z\{0}.  This  implies  that  if  A  >  1,  <p{x)  decreases  at  least 
like  |aj-2  when  x  —  x.  tan  elaborate  proof  of  a  similar  claim 
can  be  found  in  [34,  App.  B]). 

The  multiple  zeros  of  P{e 3M)  at  the  origin  [or,  equivalently, 
the  singularities  of  p(us)]  have  another  desirable  effect.  They 
will  endow  the  scaling  function  <p(x)  with  a  corresponding  order 
of  approximation  1  +  log?  A,  which  may  even  be  noninteger 


A  corresponding  (highpass)  wavelet  filter  H(e3M)  can  then  be 
determined  by  using  any  of  the  standard  design  techniques.  For 
a  complete  specification  of  the  wavelet  transform  algorithm, 
we  also  need  the  complementary  analysis  filters  G(e 3M)  and 
H  (c3"' ).  These  are  determined  by  imposing  the  perfect  recon¬ 
struction  conditions  in  the  Fourier  domain. 

Since  the  underlying  filters  are  not  necessarily  compactly 
supported,  the  best  way  to  implement  this  kind  of  wavelet  trans¬ 
form  is  to  perform  the  computation  in  the  frequency  domain.  For 
this  purpose,  we  may  use  the  fast  FFT -based  implementation2 
of  the  wavelet  transform  described  in  [35]  and  [36], 


H.  Example  1:  Fractional  Splines 

As  a  first  example  within  our  general  class  of  self-similar 
functions,  we  consider  the  first  term  (n  —  0)  of  (17):  xf  with 
a  =  log?  A  (cf.  Fig.  3).  The  Fourier  transform  of  this  one-sided 
power  function  is  given  by  the  first  term  in  (33).  The  application 
of  Theorem  3  [cf.  (35)]  leads  to  the  localization  filter 


P(e3W) 


(1  -  e~3UJ)a+1 

r(a  +  1) 


which  is  essentially  a  fractional  iterate  of  the  finite  difference 
operator  A  < — >  (1  —  e~3M).  Interestingly,  this  localization 
process  yields  the  fractional  B-splines  of  degree  a 


pl(x)  =  A“+1 


+  rfgi+ij 


Fourier 

i - - - > 


which  is  a  recently  proposed  extension  of  the  polynomial 
B-splines  for  noninteger  degrees  [34],  As  can  be  seen  in 
Fig.  4,  these  functions  provide  a  smooth  transition  between 
the  conventional  B-splines.  Each  of  them  generates  a  valid 
Riesz  basis  and  satisfies  the  partition  of  unity,  provided  that 
a  >  —(1/2)  (cf.  [34])  or,  equivalently,  A  >  l/\/2,  which  is  a 
necessary  condition  for  Theorem  2. 


2MATLAB  and  Java  sources  are  available  at  http://bigwww.epfl.ch/blu/fract- 
splinewavelets/. 
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Fig.  3.  Central  basis  functions  j>  f.r  )  =  of  the  fractional  B-splines  (see 
[37])  shown  in  Fig.  4. 


Next,  we  apply  (22)  to  determine  the  corresponding  refine¬ 
ment  filter 


G(e>w) 


1  (1  -  e-j2M)a+1 
A  (1  -  c~i^)a+1 


(1  +  e-jM)a+1 
2“ 


(23) 


which  happens  to  be  the  generalized  version  of  the  celebrated 
binomial  filter.  This  clearly  shows  that  the  fractional  B-splines 
are  valid  scaling  functions.  They  have  led  to  interesting  gen¬ 
eralizations  of  the  traditional  spline  wavelets — including  the 
Battle-Lemarie  orthogonal  splines  [35],  [36].  The  remarkable 
feature  here  is  that  these  wavelets  are  indexed  by  a  continu¬ 
ously-varying  order  parameter  (a+ 1),  as  opposed  to  an  integer, 
as  is  usually  the  case  in  wavelet  theory. 


I.  Example  2:  Harmonic  Splines 

The  second  type  of  elementary  component  in  (17)  is  of  the 
form 

^Iog2  \+j(2nir/  log  2) 

=  £+g2  A(cos(27tn  log,  x+)  +  j  sin(27rn  log,  #+)). 


These  are  essentially  modulated  versions  of  the  power  functions 
with  the  period  of  the  harmonic  component  increasing  logarith¬ 
mically.  We  can  combine  two  complex  conjuguate  components 
with  index  —  n  and  n  to  obtain  a  real-valued  central  basis  func¬ 
tion.  These  are  the  generating  functions  of  what  we  define  as 
harmonic  splines. 

Fig.  5  illustrates  the  result  of  the  localization  of  the  function 


P+(?) 


3-i(27r/log2) 


r 


3+i(27r/log2) 


r 


The  Fourier  transform  of  the  corresponding  localization  filter, 
as  given  by  Theorem  3,  is 


P(ejLJ)  = 


(1  —  e-ju^4-j(2n/  log 2) 


(1  —  e-7^)4+y(27r/  log  2) 


Fig.  4.  Fractional  B-splines  -‘r  (■>' )  with  a  >  0.  These  functions  interpolate 
the  conventional  B-splines  that  are  represented  using  a  thicker  line.  The 
noninteger  degree  splines  are  not  compactly  supported,  although  they  decrease 
like  |.r|-"-2. 


Note  that  by  construction,  p+(x)  —  8p+(x/2).  Obviously,  even 
if  its  integral  does  not  vanish  [in  fact,  f  <p(x)  dx  =  1],  this 
scaling  function  has  a  bandpass  character  that  is  rather  unusual 
for  a  scaling  function.  It  can  also  yield  wavelet  bases.  In  fact,  we 
used  the  reconstruction  part  of  a  wavelet  transform  algorithm  to 
graph  the  function  in  Fig.  5. 

IV.  From  Wavelets  to  Central  Basis  Functions 

We  will  now  show  that  one  can  also  follow  the  reverse  path 
and  uncover  the  central  basis  function(s)  that  lies  hidden  behind 
any  scaling  function  <p(x).  To  this  end,  we  need  to  assume  that 
p  is  one  sided.  This  includes  all  compactly  supported  scaling 
functions  but  excludes  those  that  are  supported  on  the  whole 
real  line. 

A.  Finding  the  Central  Basis  Function 

Theorem  4:  Let  p(x)  be  an  admissible  scaling  function  with 
corresponding  causal  refinement  filter  G{z)  —  X^>o  9kz~k 
with  <7o  0-  Then,  there  exists  a  unique  self-similar  one-sided 

(i.e.,  supported  in  [0,  +oo[)  central  basis  function  p+(x)  that 
generates  the  same  multiresolution  analysis.  The  self-similarity 
parameter  A  defined  by  (16)  is  simply  gfl . 

This  central  basis  function  p+(x)  can  be  expressed  as 

P+(x)  =  ^2  rntp(x  -  n)  (24) 

n>0 


where  the  coefficient  sequence  rn  satisfies 


To 


ffo 


Y;  ffn- 20’ k, 
k>  0 


Vn  >  0. 


(25) 


Proof:  Because  G  is  a  causal  filter,  <p  is  one-sided,  i.e., 
its  support  is  contained  in  [0,  +oo[.  We  now  dilate  the  function 


Fig.  5.  Top:  Central  basis  function  P-\-(x)  and  its  localized  version  <p(x)  (harmonic  spline).  Bottom:  Fourier  transforms  /3+(cj)  and  <p(u).  See  Section  III-I  for 
more  details. 


p ,  as  defined  by  (24),  and  apply  the  two-scale  equation  (12)  to 
show  that  it  is  self-similar. 

p{x/ 2)  =  ^2  rnip(x/2  -  n) 

n 

-  ^2  9k-2nrn<p(.x  ~  k)  [using  (12)] 

n,  k 

-  X]  9ork<p(x  -  k)  [using  (25)] 

k 

=  9op(z)- 

Note  that  we  need  not  have  technical  convergence  conditions 
for  interverting  sum  signs  here;  this  is  because  the  summation 
in  (24)  is  finite  for  fixed  x.  Thus,  we  have  span{p(a;  —  n)}  C 
span{(^(^  —  n)}. 

Conversely,  as  shown  later,  (24)  can  be  inverted  exactly  to 
yield  ■■fix)  as  a  linear  combination  of  shifted  versions  of  p+. 
This  shows  that  span{</5(/r  —  n)}  C  span{p+(a;  —  n)}.  In  fact, 
we  even  have  ip(x)  —  p+(x)  for  x  6  [0,  1[,  which  shows  that 
p. |_  is  unique  since  a  self-similar  function  is  uniquely  defined  by 
its  value  in  [1/2,  1[.  ■ 

Corollary  1:  Let  <p(x)  be  a  symmetric  compactly  supported 
scaling  function  with  corresponding  refinement  filter  G(z). 


Then,  there  exist  a  radial  basis  function  p*(x)  that  generates 
the  same  multiresolution  analysis 

p*{x)  =  p+(x)  +  P~(x)  (26) 

where  p+{x)  is  given  by  (24). 

B.  One-Sided  Localization  Filter 

In  the  present  case,  tp  is  given  a  priori ,  which  means  that 
the  localization  filter  p  is  specified  implicitly.  To  get  its  explicit 
form,  we  consider  the  matrix  form  of  (24) 


'  P+(x)  ' 

>0  ?T  ••• 

0  r0  ?’i 

o 

-  <p(x)  ' 

p+(x  -  1) 

= 

<p(x  -  1) 

(27) 

R 


where  the  infinite  matrix  R  is  Toeplitz  upper  triangular.  To  de¬ 
termine  p  such  that  <p{x)  =  [>//,> (l  pnp(x  —  n),  we  simply  need 
to  evaluate  the  inverse  of  R,  which  turns  out  to  be  Toeplitz  upper 
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0  0.5  1  1.5  2  2.5  3  3.5  4  4.5  5 


Fig.  6.  Scaling  function  ip  (dotted  line)  and  self-similar  central  basis  function 
[>  (solid  line)  corresponding  to  Daubechies  minimum  phase  wavelets  of  order 
2.  Note  that  p  and  p  coincide  over  [0,  1], 


triangular  as  well.  By  solving  this  triangular  system  of  equa¬ 
tions,  we  find  that 


k 

^  [ knPk—n' 
n= 1 


(28) 


Note  that  this  localization  is  not  necessarily  the  same  as  the 
one  described  in  Theorem  3,  which  may  yield  another  scaling 
function  in  V(l . 


C.  Example  1:  Daubechies 

Daubechies  scaling  functions  [38]  are  compactly  supported 
and  can  thus  be  expressed  using  a  unique  self-similar  one-sided 
central  basis  function.  To  illustrate  our  theory,  we  have  eval¬ 
uated  the  sequence  rn  using  (25)  for  a  Daubechies  filter  of 
length  4,  which  corresponds  to  a  wavelet  of  order  2.  The  cor¬ 
responding  functions  are  plotted  in  Fig.  6.  Note  that  p+,  which 
is  clearly  fractal,  also  corresponds  to  the  self-similar  extrapola¬ 
tion  of  </2[i/2,  i] ,  i.e.,  the  restriction  of  <p{x)  to  x  6  [1/2,  1],  with 
A  =  4/(1  +  v%)  ~  1-46. 


One  can  check  that  the  results  for  n  —  1  are  compatible  with 
(7);  therefore,  we  have  gone  full  circle. 

E.  Time-Domain  Parametrization  of  Scaling  Functions 

Theorem  4  allows  us  to  state  the  utmost  result  of  this  paper, 
i.e.,  a  full  parametrization  of  one-sided  dyadic  scaling  functions. 

Corollary  2:  Let  ip  be  an  admissible  one-sided  or  compactly 
supported  scaling  function.  Then,  there  exist  a  real  number  a  > 
—0.5,  a  complex  sequence  {7 n}nei  £  and  a  real  sequence 
{Pnjne N  such  that 

p(x)  -  ^2  -  k)++J('2l7r/los2) ■  (31) 

fc€N,  lez 

Conversely,  any  function  described  by  (31)  satisfies  a 
two-scale  difference  equation  but  is  not  necessarily  compactly 
supported. 

As  far  as  we  know,  this  is  the  first  time  that  an  explicit  time 
domain  formula  is  given  for  scaling  functions.  Equation  (31)  is 
interesting  because  it  explains  the  fractal  nature  of  scaling  func¬ 
tions  and  wavelets.  The  expansion  can  also  be  interpreted  as  an 
infinite  sum  of  harmonic  splines.  This  points  to  a  new  funda¬ 
mental  aspect  of  splines  as  elementary  constituents  of  scaling 
functions. 

Note  that  the  exponent  a  in  (31)  has  a  very  special  role. 
Whenever  the  sum  over  7 ;  is  finite,  it  is  tightly  linked  to  four 
important  wavelet  properties: 

i)  the  Holder  regularity,  which  is  (  v; 

ii)  the  Solobev  regularity,  which  is  a  +  (1/2); 

iii)  the  order  of  approximation,  which  is  a  +  1; 

iv)  the  polynomial  degree  reconstruction,  which  is  [a]. 

We  remark  that  except  for  the  case  of  splines,  there  is  a  whole 
range  of  such  wavelets  that  has  not  yet  been  explored. 

We  should  mention,  however,  that  the  meaning  of 
a  —  —  log2  <7o  is  not  as  transparent  for  the  traditional 

wavelets  such  as  Daubechies  whose  harmonic  decomposition 
involves  an  infinite  number  of  terms.  This  leaves  us  with  the 
following  interrogation:  The  Holder  regularity  of  the  minimum 
phase  Daubechies  wavelets  turns  out  to  be  a  —  —  log2  go  (cf. 
[39]— [41]).  Is  this  fact  coincidental  or  not? 


D.  Example  2:  Polynomial  Splines 

The  same  kind  of  exercise  can  be  performed  for  the  causal 
B-splines  (15)  to  recover  the  one-sided  power  functions  xf. 
Here,  instead  of  using  (25),  we  will  determine  the  sequence  r 
as  the  inverse  of  the  localization  operator  p  to  illustrate  the  link 
in  Section  II-A. 

The  ^-transform  of  the  localization  operator  for  the  B-spline 
of  degree  n  is 

Pn(z)=^j(l-z-l)n+l.  (29) 

Formally,  11"  (z)  is  obtained  by  taking  the  convolution  inverse  of 
Pn(z).  Using  the  Taylor  development  of  (1  —  z-1)-"-1  (which 
is  given,  e.g.,  in  [28]),  this  yields 


F.  Borderline  Cases 

We  have  established  a  full  equivalence  between  self-similar 
central  functions  and  scaling  functions  only  for  the  one-sided 
case.  It  is  likely  that  a  similar  equivalence  holds  for  most  scaling 
functions  that  are  not  one-sided  as  well,  but  we  also  have  evi¬ 
dence  that  this  will  not  always  be  the  case. 

An  interesting  counterexample  is  provided  by  the  sym¬ 
metric  fractional  splines  whose  generating  function  is 
(>,,  (x)  =  x2n  log  x  | ,  where  n  is  some  integer.  We  have  shown 
recently  that  this  classical  radial  basis  function  could  be  local¬ 
ized  to  yield  the  symmetric  fractional  B-spline  fl2n,  which  is 
a  perfectly  valid  scaling  function  [34].  Yet,  the  corresponding 
(>+  (x)  is  not  strictly  self-similar;  it  is  only  self-similar  up  to 
a  polynomial  term:  p*(x/ 2)  =  {1 / An)(p*(x)  —  a:2" log 2). 
However,  the  span  of  these  radial  basis  functions  is  still  a 
multiresolution  space  because  x2n  belongs  to  it. 
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V.  Conclusion 

We  have  presented  new  results  that  make  the  connection  be¬ 
tween  radial  basis  functions  and  wavelets  very  explicit. 

The  fact  that  we  can  move  from  radial  basis  functions  to 
wavelets  enables  us  to  control  some  of  their  key  mathemat¬ 
ical  properties:  regularity  and  order  of  approximation.  It  also 
yields  wavelets  that  have  an  explicit  analytical  form,  the  frac¬ 
tional  splines  being  a  notable  example. 

The  existence  of  a  link  in  the  reverse  direction — from 
wavelets  to  central  basis  functions — is  especially  interesting 
conceptually.  It  leads  to  an  alternative  interpretation  of  mul¬ 
tiresolution.  Basis  functions  are  simply  removed  (resp.  added) 
instead  of  being  dilated  (resp.  contracted),  as  is  usually  the 
case.  This  opens  up  the  door  to  many  possible  extensions,  such 
as  wavelets  on  nonuniform  grids. 

This  study  has  also  brought  out  the  fundamental  role  of 
splines  and  their  harmonic  extensions  that  are  the  elementary 
constituents  of  wavelets.  The  decomposition  in  terms  of  har¬ 
monic  splines  leads  to  an  explicit  formula  for  scaling  functions 
and  wavelets  in  the  continuous-time  domain.  It  also  provides  a 
clear  understanding  of  their  fractal  properties. 


Appendix  A 

Fourier  Transform  of  a  One-Sided  Self-Similar 
Central  Basis  Function 

It  turns  out  that  the  Fourier  transform  of  the  tempered  distri¬ 
bution  +  (x)  —  x'p  is  known  in  the  general  case  where  t£C. 
It  is  given  by  (see  [42]) 

u(to)  —  if  x  N 

V  ’  (MA+1  * 

u(u)  =  ,.n',  +  7t jn6^(w)  if n  £  N. 

Thus,  in  the  sense  of  distributions,  we  have 


P^)=^lnT. 


r  a  +  1  +  j 


2nir 

iog2 


nez 


(jo;)a+1+f  (2,!7r/  los  2 


£(w)  =  Z)  7’!  T7 


r  a  +  1  +  j 


2mr 

bg2 


nez 


(jw)a+1+-,’(2mr/  loS  2) 


if  a  —  log,  A  6 


if  a  =  log,  A  f 

■707 xjaS(a\u) 


(32) 


Appendix  B 
Proof  of  Theorem  3 

Proof:  Here,  we  assume  that  p  can  be  expressed  as  in  (17) 
with  a  finite  sum.  From  Appendix  A,  we  also  have  its  Fourier 
transform  in  the  sense  of  distributions: 

_  A  ,  .  .  2mr  A 

r  ( i  +  log,  a  +  j  - — -  ) 

Y  x  _  _V _ log  2  J  m 

Pv0)  —  2—t  \l+log2  A+/(2mr/log2) 

finite  n  V  / 

if  A  is  not  an  integer  power  of  two.  In  the  case  where  A  =  2k  for 
some  k  £  N,  (u>)  has  to  be  added  to  this  expression. 


We  now  use  (20)  to  derive  the  filter  contribution  associated 
with  a  single  component  of  the  form  xf,  making  the  complex 
change  of  variable3  u  —  (1  —  z_1)x 


cae-^  1>dx  = 


1-z 


-l 


du 


1-z 


-l 


(1  -*-1)a+1 
T(1  +  a) 

(1  -z-1)^1' 


(34) 


By  summing  over  all  components,  we  get  the  explicit  form  of 
the  inverse  of  the  localization  filter 


F(e^)-1 


r(i  +  iog,A  +  j^) 

v  7  _ A _ log>2/ _  (35) 

n  —  e_7^)l  +  log2  A+j(2)!7 t/  log  2) 
finite  n  ^  ' 


Because  1  —  e~JW  —  ju>  +  0{ur),  it  is  clear  that  P{eJW)~l  is 
equivalent  to  p(u>)  in  the  neighborhood  of  to  —  0  [except  for  the 
potential  S^(ui)  correction  in  the  few  cases  where  A  =  2k]. 
Clearly,  the  localization  filter  l’ic:'"  )  will  be  well  defined  and 
bounded  provided  that  P(e-,w)_1  0  over  [0,  27t[.  Under  those 

conditions,  we  can  state  that  <p{u>)  —  P(e-,M)p(oj)  is  bounded 
over  FS  due  to  the  fact  that  p(u>)  is  bounded  away  from  u>  —  0. 
The  same  argument  holds  for  A  =  2k  because  then,  P(cJ^)  — 
0(u>k+1)  and  because  wfc+1 6^  (w)  =  0.  Now,  we  claim  that  (p 
is  in  L2(R).  This  is  because 


1 

(ja;)1+lo82  A+y(2)!7r/  log  2) 


e)!7T2/log2sign(a;) 
1^11+7^(1082  A) 
Constant 
|a;|l+^(log2  A) 


which  implies  that  \  0{lo)\2  <  Constant  X  |w|  2  27^(lo82  A)  away 
from  to  —  0.  Thus,  \(p\ 2  is  integrable  over  ]— oo,  -i]  U  [1  ,  oo[, 
provided  that  |A|  >  1/V2.  Since  \(p\ 2  is  also  bounded,  it  is 
obviously  integrable  over  [—1,  1].  By  Plancherel  identity,  we 
conclude  that  ip  £  L2(R). 

If  p  is  symmetric,  we  decompose  it  as  p+(x)  +  p+(—x)  and 
proceed  similarly  to  show  that  <p{u>)  —  P(e3M)p(oj)  is  bounded 
and  square  integrable,  which  proves  that  ip  is  square  integrable. 
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